/* es-replication.do --- 
 * 
 * Filename: es-replication.do
 * Description: (Hopefully) replicates the analyses in my ES piece
 * Author: Kai Arzheimer
 */

/* Commentary: 
 * Writing the regression tables requires outreg2
 * Install it with ssc install outreg2
 * 
 */

/* Change log:
 * 
 * 
 */

/* Code: */

* Set mem, load data, generate dependend variable

set mem 90m
use pi-west-77-02,clear
ren pi piorig
recode piorig 1/4=1,gen(pi)


* Create (Dummy)-Variables
gen kath=0
replace kath=1 if konf==2
replace kath=. if konf==.

/*Min Kirchgang und Zeit =0*/
gen nullkirch=kirchgang -1
gen nullzeit=zeitp-206

/*Interaktionen*/
gen kirchXzeit=nullkirch*nullzeit
gen arbXzeit=nullzeit*arbeiter
gen kathXzeit=nullzeit*kath

gen kathXarb=arbeiter*kath
gen kathXarbXzeit=kathXarb*nullzeit
gen kathXkirch=kath*nullkirch
gen kathXkirchXzeit=kath*nullkirch*nullzeit


* Graph No1
preserve
* Non-missing cases
gen piok= pi~=.
bysort zeitp: egen pioksum=sum(piok)
replace pi=pi*100
collapse pi pioksum,by(zeitp)
summ pioksum


gen rssize=pioksum/985
/*divide by size of average sample*/
gen sqrrssize=sqrt(rssize)
/*Take square root*/

twoway (scatter pi zeitp [aw=rssize], mcolor(black) msymbol(oh) mcolor(gs8)) (lowess pi zeitp,clpattern(solid) ), plotregion(style(none)) legend(order(1 "monthly samples" 2 "lowess smoother (bw=0.8)" ) rows(1)) xtitle("") ytitle("Percent Identifiers")  yscale(r(54 84) extend) ylabel(54(10)84) note("symbols scaled proportional to sample size")

graph save fig1,replace
restore

* Table 2: Models
	/*Die Zeit-Variablen f�r die Schnittpunkte*/

gen postunification=0
replace postunification=1 if jahr>1990

gen ab92=0
replace ab92=1 if jahr >1991

gen ab82=0
replace ab82=1 if jahr >1981
	
gen postXnullzeit=postunification*nullzeit
gen ab92Xnullzeit=ab92*nullzeit
gen ab82Xnullzeit=ab82*nullzeit

logit pi nullzeit ,cluster(zeitp)
predict ppur

	
outreg2 using tabelle2.txt,coefastr se replace 

logit pi nullzeit postunification postXnullzeit ,cluster(zeitp)
predict p91
	
	
outreg2 using tabelle2.txt,coefastr se append 

logit pi nullzeit ab92 ab92Xnullzeit,cluster(zeitp)
predict p92
	
outreg2 using tabelle2.txt,coefastr se append
	

logit pi nullzeit ab82 ab82Xnullzeit,cluster(zeitp)
predict p82


outreg2 using tabelle2.txt,coefastr se append 

* Figure 2
	/*Und hier noch mal die Grafik*/
foreach var of varlist ppur p91 p92 p82 {
	replace `var'=`var'*100
	}

preserve

collapse ppur p91 p92 p82,by(zeitp)

twoway (line ppur p91 p92 p82 zeitp,clpat(solid longdash shortdash dot) clcolor(black gs7 gs7 gs7)), plotregion(style(none)) legend(order(1 "no break" 2 "break '91" 3 "break '92" 4 "break '82") size(small) subtitle("alternative trend specifications",size(small) ) rows(1)) xtitle("")   yscale(r(54 84) extend) ylabel(54(10)84)  ytitle("Percent Identifiers")

graph save fig2,replace

restore


* Table 3 contra Dalton
gen abiXzeit=abi*nullzeit
logit pi abi nullzeit abiXzeit ,cluster(zeitp)
outreg2 using tabelle3.txt,coefastr se replace

* Table 4: the monster

logit pi nullzeit arbeiter kath nullkirch  kathXarb  kathXkirch arbXzeit kathXzeit kirchXzeit  kathXarbXzeit kathXkirchXzeit if piori<3,cluster(zeitp)
	
outreg2 using tabelle4.txt,coefastr se replace

* Figure 3
gen simzeit=.

/*Worker, catholic, church attendance 3*/
gen arbkath3=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace arbkath3=100*(exp(_b[_cons]+_b[arbeiter]+_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+_b[kathXarb]+3*_b[kathXkirch]+`j'*3*_b[kathXkirchXzeit]+`j'*_b[kathXarbXzeit])/(1+exp(_b[_cons]+_b[arbeiter]+_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+_b[kathXarb]+3*_b[kathXkirch]+`j'*3*_b[kathXkirchXzeit]+`j'*_b[kathXarbXzeit]))) in `i'
	}


/*Worker, non-catholic, church attendance 3*/

gen arbnonkath3=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace arbnonkath3=100*(exp(_b[_cons]+_b[arbeiter]+0*_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+0*3*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+_b[arbeiter]+0*_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+0*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+0*_b[kathXarbXzeit]))) in `i'
	}

	
/*Worker, non catholic, no church attendance*/

gen arbnonkath0=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace arbnonkath0=100*(exp(_b[_cons]+_b[arbeiter]+0*_b[kath]+0*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+`j'*0*_b[kathXzeit]+`j'*0*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+0*`j'*_b[kathXarbXzeit])/(1+exp(_b[_cons]+_b[arbeiter]+0*_b[kath]+0*_b[nullkirch]+`j'*_b[nullzeit]+`j'*_b[arbXzeit]+0*`j'*_b[kathXzeit]+`j'*0*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+0*`j'*_b[kathXarbXzeit]))) in `i'
	}

/*Non-workers, catholic, max church attendance*/

gen nonarbkath5=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace nonarbkath5=100*(exp(_b[_cons]+0*_b[arbeiter]+_b[kath]+5*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*5*_b[kirchXzeit]+0*_b[kathXarb]+5*_b[kathXkirch]+`j'*5*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+0*_b[arbeiter]+_b[kath]+5*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*5*_b[kirchXzeit]+0*_b[kathXarb]+5*_b[kathXkirch]+`j'*5*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit]))) in `i'
	}

/*Non-workers, non-catholic, max church attendance*/

gen nonarbnonkath5=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace nonarbnonkath5=100*(exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+5*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+0*_b[kathXzeit]+`j'*5*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+5*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+0*_b[kathXzeit]+`j'*5*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit]))) in `i'
	}

/*Non-workers, non-catholic, church attendance 3*/

gen nonarbnonkath3=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace nonarbnonkath3=100*(exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+0*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+0*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit]))) in `i'
	}

/*Non-workers, catholic, church attendance 3*/

gen nonarbkath3=.

forvalues j = 0(1)309 {
	local i= `j'+1
	replace simzeit=`j' in `i'
	replace nonarbkath3=100*(exp(_b[_cons]+0*_b[arbeiter]+_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+3*_b[kathXkirch]+`j'*3*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+0*_b[arbeiter]+_b[kath]+3*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*_b[kathXzeit]+`j'*3*_b[kirchXzeit]+0*_b[kathXarb]+3*_b[kathXkirch]+`j'*3*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit]))) in `i'
	}

	/*Reference: Non-workers, non-catholics, non-church-goers*/

gen nonarbnonkath0=.

forvalues j = 0(1)309 {
	local i= `j'+1

	replace nonarbnonkath0=100*(exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+0*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*0*_b[kathXzeit]+`j'*0*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit])/(1+exp(_b[_cons]+0*_b[arbeiter]+0*_b[kath]+0*_b[nullkirch]+`j'*_b[nullzeit]+`j'*0*_b[arbXzeit]+`j'*0*_b[kathXzeit]+`j'*0*_b[kirchXzeit]+0*_b[kathXarb]+0*_b[kathXkirch]+`j'*0*_b[kathXkirchXzeit]+`j'*0*_b[kathXarbXzeit]))) in `i'
	}

	
replace simzeit=simzeit+206
format %tm simzeit
	
twoway (line arbkath3 arbnonkath0 nonarbkath5 nonarbnonkath5 nonarbnonkath3 nonarbkath3 nonarbnonkath0 simzeit,clpat(dash dot solid solid dash dash dot) clcolor(gs7 black black black gs7 gs7 black)), plotregion(style(none)) legend(order(2 "never" 1 "occasionaly" 3 " every Sunday" ) size(small) subtitle("church attendance",size(small) ) rows(1)) xtitle("")   yscale(r(54 84) extend) ylabel(54(10)84)  ytitle("Percent Identifiers")  text(73  190 "-catholic /+worker",place(e) size(vsmall)) text(69.5  195 "-catholic /-worker (ref)",place(e) size(vsmall)) text(84  216 "+catholic /-worker",place(e) size(vsmall)) text(69.5  195 "-catholic /-worker (ref)",place(e) size(vsmall)) text(68.5  460 "-catholic /-worker",place(e) size(vsmall)) text(61.8  460 "+catholic /+worker",place(e) size(vsmall)) text(75.5  300 "+catholic /-worker",place(e) size(vsmall)) text(76.4  175 "-catholic /-worker",place(e) size(vsmall)) 

graph save fig3,replace


/* es-replication.do ends here */
